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Crucially important for application of type-II superconductor films is the stability of the vortex matter - 
magnetic flux lines penetrating the material. If some vortices get detached from pinning centres, the energy 
dissipated by their motion will facilitate further depinning, and may trigger a massive electromagnetic 
breakdown. Up to now, the time-resolved behaviour of these ultra-fast events was essentially unknown. We 
report numerical simulation results revealing the detailed dynamics during breakdown as within 
nanoseconds it develops branching structures in the electromagnetic fields and temperature, with striking 
resemblance of atmospheric lightning. During a dendritic avalanche the superconductor is locally heated 
above its critical temperature, while electrical fields rise to several kV/m as the front propagates at instant 
speeds near up to 100 km/s. The numerical approach provides an efficient framework for understanding the 
ultra-fast coupled non-local dynamics of electromagnetic fields and dissipation in superconductor films. 

The smooth gradual penetration of magnetic flux in type-II superconductor films experiencing an increasing 
perpendicular magnetic field can be interrupted by large amounts of flux rushing in almost instantaneously. 
Left to be observed after such an event is the intriguing dendritic structure frozen in the flux distribution, see 
Ref. 1 for a review. This dramatic avalanche phenomenon has been observed in films of many materials, such as 
Nb 2 ' 3 , YBa 2 Cu 3 O x 4 - 7 , Pb 8 , MgB 2 9 " 15 , Nb 3 Sn 16 , YNi 2 B 2 C 17 , NbN 1819 , and very recently also in a-MoGe 20 . It is today 
widely accepted that such events originate from a thermoelectromagnetic instability 21 , where an initial fluctuation 
sets some flux in motion and the dissipative displacement increases the local temperature. This reduces the flux 
pinning, and facilitates even further motion of flux in the same region. Thus, a positive feedback loop is created, 
which may lead to a sudden catastrophic breakdown in the physical properties, and hence, also the potential 
practical use of the superconductor. 

These events are extreme also because they occur unpredictably and develop too fast to be followed dynam- 
ically by any experimental method available today. With conventional magnetometry one observes only a step in 
the magnetic moment due to the abrupt redistribution of flux and induced currents 22 ' 23 . More information is 
obtained from magneto -optical imaging (MOI), where the spatial distribution of magnetic flux before and after 
the breakdown is visualized using a Faraday- active sensor mounted on the sample 24 . However, results providing 
insight into how the breakdown evolves in time are extremely scarce. Only by using a femtosecond pulsed laser to 
actively trigger an event it was possible to synchronize the image recording and to capture the flux distribution at 
an intermediate stage 5 " 7 . 

From those experiments it was found that the flux front advances in a nonballistic way with a velocity decaying 
in time. In the initial stage the front can advance at an astonishing speed of 100 km/s. This dramatic dynamics 
poses a long list of unresolved questions, in particular, (i) How do the electric field, temperature, and flow of 
shielding currents evolve during the breakdown? (ii) In a typical event many branches are formed. Do they grow 
simultaneously, or develop one by one from the same root? (Hi) Where is the temperature maximum, at the tip or 
at the root of the dendritic structure, and does it exceed the critical temperature T c of the superconductor? And 
(iv) Does there exist any physical activity that could serve as a warning that breakdown is near? 

We have addressed these questions by performing numerical simulations of the thermo- electromagnetic 
behavior of a superconducting thin disk. The modelling is done using classical electromagnetism with a highly 
nonlinear material law 25 coupled to an equation describing the production and transfer of heat. In films, com- 
pared to bulk superconductors, the theoretical treatment is far more challenging due to the non-local electro- 
dynamics and the presence of a substrate at constant temperature T 0 in close thermal contact with the sample 26 . 

Linear stability analysis starting from a smooth critical- state flux distribution has provided the conditions for 
onset of the thermomagnetic instability. The most unstable mode has a finite wave-number, which means that 
the instability will nucleate by fingering 27 " 31 . As the instability develops, numerical analysis of the non-linear and 
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Figure 1 | Sample configuration. A thin superconducting disk on a 
substrate exposed to a gradually increasing perpendicular magnetic field, 
H a . The flux density, B z , is advancing from the edge along with a 
distribution of induced shielding current, /, and electrical field, E. 

non-local dynamics of magnetic flux and temperature have revealed 
that the instabilities develop into complex branching structures 28,32 . 

Here we focus on the temporal evolution of the instability - a 
complete scenario beyond experimental accessibility. Our results 
cover the nucleation stage of the breakdown, the branching stage 
when the dendritic pattern emerges, and finally the stagnation of 
flux motion and decay of the released heat. We explain the physical 
mechanisms involved at every stage of the event, and provide ana- 
lytical expressions for all the relevant time scales. 

Results 

Flux penetration. The simulations were carried out for a thin disk 
geometry (Fig. 1) using material parameters corresponding to 
superconducting MgB 2 . The Methods section gives a description of 
the numerical approach. A magnetic field H a is applied transverse to 
the sample plane, and as it gradually increases from zero it drives the 
penetration of magnetic flux into the disk. In the early stage of the 
field ramp, the flux enters evenly around the edge, and advances to 
increasing depth without any sign of intermittent behaviour. In the 
penetrated region a critical state is formed and characterized by a 
sheet current / and flux density B z in full agreement with the Bean 
model for a thin circular disk 33 " 35 . 

When the applied field reaches ^ 0 H th = 5.3 mT the first abrupt 
event is nucleated, and magnetic flux rushes in from the edge; see 
Supplementary Movie 1. A complex branching structure is created as 
the flux invades deep into the Meissner state region, which previously 
was free from flux, see Fig. 2a. As H a continues increasing, only the 
gradual flux penetration resumes while the dendritic structure 
remains frozen. Then, at the field of 6.2 mT, another similar event 
takes place in a different part of the sample, and soon thereafter yet 
another one strikes. In this way the superconductor experiences a 
sequence of dramatic events at unpredictable intervals and locations, 
and where each breakdown follows an intriguing path in a macro- 
scopically uniform medium. Since this phenomenon is of electro- 
dynamic nature, it is interesting to recognize the many aspects that 
are similar to atmospheric lightening. One exception, however, is 
that in the superconductor the magnetic footprints are preserved, 
which has allowed experimental visualization of the phenomenon. 
Figure 2b shows an MOI picture of the flux distribution in a super- 
conducting MgB 2 film at T = 5 K where the magnetic field had been 
increased from zero to id 0 H a = 3.8 mT. The experimental image 
reveals that the flux avalanches have a morphology quite similar to 
the numerical results, and also that the events have a clear tendency 
to avoid spatial overlap, as in the simulations. 

Evidently, in a circular disk there is no geometrically preferred 
location for the first breakdown to nucleate. Instead, the place of 
the first event is selected by a combination of the quenched disorder 




Figure 2 | Flux density after a few breakdown events, (a) Simulated 
distribution of B z in a superconducting disk after five flux avalanches 
occurred in the sequence indicated by the numbers as the applied field was 
ramped up from zero to fi 0 H a = 8.5 mT. (See Supplementary Movie 1.) (b) 
Magnetooptical image of the flux density in a superconducting MgB 2 film 
cooled to 6 K and then exposed to an applied field of 3.8 mT. 

and the random fluctuations. However, after the first event has 
occurred, the symmetry is broken because the flux avalanche 
destroys locally the critical -state, and thus reduces the shielding cur- 
rent near the edge where the flux enters. As a result, typical for films 
placed in a perpendicular magnetic field, the largest edge field is 
found in sectors between the roots of previous events. Hence, there 
is an anti-correlation in space and time between avalanches. The 
nucleation sites therefore tend to distribute evenly around the cir- 
cumference, as one finds also experimentally (Fig. 2b). 

Figure 3 shows the total magnetic flux O = JcPr B z (r) integrated over 
the disk area plotted as function of the increasing applied field. Initially, 
the 0>(H a ) grows gradually until the threshold field of 5.3 mT is reached. 
Here <I> makes a jump corresponding to entrance of 2.6 millions flux 
quanta, and forms the dendritic structure labelled (1) in the image. The 
subsequent jumps marked (2) -(5) involve a similar amount of flux. The 
flux entering during the intervals between subsequent jumps is also 
substantial, but smaller than the flux involved in each event. 

Electric field and temperature. Far less accessible by experimental 
methods is the behaviour of the electric field E and temperature T in 
the superconductor, and we report here only results obtained from 
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Figure 3 | Magnetic flux in the sample versus field. The graph displays the 
total magnetic flux, Q, inside the sample as the applied field is increased in 
the simulations. The avalanche events are seen as a sequence of vertical 
steps in the curve. The amount of flux is measured in millions of flux 
quanta O 0 = hl2e. The five images show the flux distribution right after 
each event has occurred in the sequence (l)-(5). 



the simulations. Shown in Fig. 4a and b are the maximum values £ max 
and T max in the disk area recorded as the magnetic field is increasing. 

At small H a the £-field is fairly small, typically in the range of 0.1- 
1 V/m. Due to the disorder introduced in the model, the electric field 
fluctuates and the location of £ max is not fixed, though it is always 
found near the edge. The overall trend is that £ max steadily increases 
with H a , a behaviour which is in agreement with the Bean model 
where EccH a H a at small magnetic fields 25 . Yet, our £ max is consid- 
erably larger than the maximum value in a uniform disk, which is 
E = fi 0 H a R/2 = 0.01 V/m in the Bean model at full flux penetration 35 . 
The fairly strong fluctuations seen in the £-field up to /n 0 H a = 5.3 mT 
suggest that the flux penetration proceeds in small bursts. Such 
micro-avalanches, which locally enhance the £-field, then become 
potential nucleation sites for the largescale events. Indeed, previous 
MOI studies have reported observations fully consistent with this 
picture 36,37 . 

At fioH a = 5.3 mT the first big event strikes, leading to a spike in 
-Emax peaking at nearly 10 kV/m. Immediately after, the £-field drops 
by five orders of magnitude. Then, over a substantial ramping inter- 
val the £ max stays far below the value it had before the event. This 
means that micro -avalanches are strongly suppressed for some per- 
iod. As H a continues increasing the maximum E- field gradually 
returns to the pre-event level, and the next breakdown occurs 
1 mT after the previous one. The graph (Fig. 4a) shows that the same 
characteristic behavior is found before and after each big event. 

The maximum temperature over the disk area (Fig. 4b) is initially 
equal to the substrate temperature, and the superconductor is ther- 
mally uniform. As the magnetic field increases, precursors of a big 
event begin to appear as noise also in T max . From jd 0 H a ~ 4 mT the 
fluctuations increase in amplitude, and at 5.3 mT the temperature 
shoots up to T max = 1.3T C . Thus, during the breakdown supercon- 
ductivity is lost in part of the sample. As for the electric field, right 
after the event the fluctuations in T max are suppressed, before they 
gradually reappear as the ramping of H a continues. Similar to £ max , 
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Figure 4 | Maximum electric field and temperature, (a) Maximum local 
electric field during the magnetic field ramp. The E-field is shown in a 
logarithmic scale, (b) Maximum local temperature in the disk normalized 
to the critical temperature T c of the superconductor. The base line 
represents the temperature of the substrate, T 0 = 0.2 T c , and the distinct 
spikes correspond to the five events shown in Fig. 3. 
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is nearly the same. 
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Time evolution. Consider now the detailed dynamics of one 
breakdown, and we choose to zoom in on the event taking place at 
/n 0 H a = 5.3 mT. Shown in Fig. 5 rows (a)-(d) are five instantaneous 
distributions of the magnetic flux density B z , the stream line pattern 
of the flow of sheet current /, the temperature T, and the electric field 
£, respectively. The snapshots show the states at t = 1,5, 22, 52 and 
86 ns after nucleation of the instability. The time evolution of B z and 
T can be seen in Supplementary Movie 2. 

The invading flux, row (a), develops from being at t = 1 ns only a 
small protrusion on top of the critical state, to become a large dend- 
ritic structure after 86 ns. Already at t = 5 ns a branching process has 
started, and at t = 22 ns the second generation branches are clearly 
visible. Between t = 52 and 86 ns the structure develops less fast. The 
final flux distribution looks quite similar to those reported from 
many MOI experiments 2 " 19 . Interestingly, it is essentially impossible 
to tell from the distribution at t = 22 ns whether it is the final or an 
intermediate state of the flux avalanche. Indeed, considerable self- 
similarity in the evolving flux pattern has been reported also from 
ultra-fast MOI experiments 5 " 7 . 

In the image at t = 52 ns, the location marked (A) shows an 
example of one type of branching, namely where a side branch starts 
growing from an already existing main branch. Another type of 
branching can be seen at location (B), where the tip of a growing 
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branch suddenly splits almost symmetrically, and both new branches 
grow simultaneously ending up equally long. The simulations show 
that the latter type is the less frequent branching process. From the 
row of B z images one can also see that the branches stop developing at 
different times. The main branches grow throughout the whole event 
and become the longest, not because of an extreme growth rate but 
because they develop over a long time span. The flux density is quite 
uniform along the branches, and has values that can exceed /HoH a by a 



factor of two or even more. Note also that along the central part of 
some branches the flux density is strongly reduced. Both these char- 
acteristics reproduce features reported from previous experiments 2,9,10 . 

Row (b) shows how the current flow pattern evolves from an initial 
criticalstate flow, where the stream line density is constant in the 
flux-penetrated region, to a complex flow pattern developing as 
the avalanche proceeds. A front of high current density follows the 
expanding perimeter of the growing flux structure, leaving behind 
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Figure 6 | Velocity and dissipation, (a) Length, L, of the longest branch in 
a flux avalanche, and its propagation velocity, v p as functions of time, 
(b) Transient behavior of the heat production rate P in the disk and the 
integrated heat release Q. 

regions of sub -critical current density. Reduced current density 
inside the frozen dendritic structure has also been seen experiment- 
ally 13,38 . Note that in addition to the currents flowing along the 
branches, there is also a component of uniform current flowing 
through the entire dendritic structure. This transverse current acts 
as a driving force for the propagation of the flux avalanche. 

Row (c) displays the elevated temperature T, which has the same 
overall dendritic signature as the flux density. Already after 1 ns a 
heated region extends from the edge to the flux front. In the next 
three panels, the branching structures develop rapidly, quite similar 
to the pattern of the moving flux. The avalanche stops propagating 
when heat removal to the substrate eventually balances the heat 
production. This is seen in the final panel, where the maximum 
temperature has already decayed to 0.7 T c . Shortly after that no ther- 
mal traces are left from this dramatic event where superconductivity 
was lost for a very short period of time. 

At all stages, the T-maps give a very good view of the dynamics, 
since only the moving branches are producing heat and excess heat is 
rapidly removed due to the good thermal contact with the substrate. 
Even after t = 22 ns, some of the minor branches are relatively cold 
(dark red), indicating that here the heat generation is reduced. The 
maximum temperature is all the time found in the trunk of the den- 
dritic structure, and most of the time it exceeds T c . The tip of the 
branches are heated significantly, but the temperature is less than T c . 

Row (d) shows the distributions of the electric field, E, where the 
magnitude is displayed on a logarithmic scale. Evidently, there is a 
strong correlation between the regions with high £-field, and the 
regions of large temperature rise. One finds here E-fields up to several 
kV/m, which is consistent with p 0 j c (T 0 ) = 6.7 kV/m, as expected 
when the sample enters the Ohmic regime. This therefore suggests 
that the high electric fields are created while the supercurrents with 



critical current density j c for a short time flow locally with Ohmic 
resistivity. 

Velocity and dissipation. From the simulations one can extract also 
the rate of propagation of the dendritic structures. Let L denote the 
distance from the sample edge to the tip of the flux branch penetrated 
deepest into the disk. For the event shown in Fig. 5 one obtains the 
curve for L(t) displayed in Fig. 6a. The L(0) = 0.12.R is the overall flux 
penetration depth prior to the event. After t = 75 ns the motion is 
over, ending with L = 0.78R, which gives an average velocity of 
18 km/s. The plot includes also the instantaneous propagation velo- 
city, v p (t) = dL/dt, of the avalanche. Right after onset the speed is 
close to 90 km/s, and then decaying during 50 ns to values near 
10 km/s, a behaviour in very good quantitative agreement with 
experimental results 5 ' 6 . Note that the initial velocity exceeds the 
speed of sound, and the propagation at this stage is driven by the 
nonlocal electrodynamics and adiabatic heating rather than by 
thermal diffusion; see the Discussion section for quantitative esti- 
mates of characteristic times. 

Presented in Fig. 6b is the transient behaviour of the heat dissipa- 
tion rate, P(t) = J disk d 2 r JE. It rises to a maximum of nearly 1 W after 
approximately 20 ns. The subsequent decay is essentially a con- 
sequence of the active area of the avalanche beginning to shrink. 
This is fully consistent with the £-field maps of Fig. 5 where the 
dendritic structure at t = 22 ns is the largest in size and intensity. 
Included in the plot is also a graph of the released heat, Q(t) = 
J 0 dt' P(t'), which during the event amounts to a total of 40 nj. 

High time-resolution plots of T max and £ max during the break- 
down event are shown in Fig. 7. The transient behavior of T max is 
similar to that of the integral dissipation rate, P(t). The main differ- 
ences are that the peak temperature occurs earlier, only 10 ns 
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Figure 7 | Transient temperature and electrical field. Global maximum 
temperature (a) and E-field (b) in the disk as functions of time during the 
breakdown event. The points marked on the T-curve correspond to the five 
temperature distribution shown in Fig. 5. Note also that both these graphs 
are zooming in on the first spike in the plots of Fig. 4. 
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after instability onset, and that the temperature pulse decays slower, 
which is due to the thermal inertia of the system. Note from the graph 
that just prior to onset of the event, the T max increases by 1-2 K 
during a 2 ns period. This is consistent with our suggestion that 
micro -avalanches form the triggering mechanism of the breakdown. 
This is further supported by the behavior of the £ max in the graph 
below. There one sees a clear signature of micro -avalanches taking 
place during the last couple of nanoseconds before £ max abruptly 
jumps by four orders of magnitude. Interestingly, the £-field also 
reveals that a very short-lived event takes place approximately 
10 ns after the main event terminated. 

Discussion 

The present simulations have revealed several important time scales 
characterizing the nucleation and subsequent evolution of the thermo- 
electromagnetic breakdown in type- II superconducting films. Most 
notably, we find that the rise time of such events, described by how 
fast the electric field rises to its maximum (Fig. 7b), is extremely short: 
of the order of 1 ns. The total duration of an event is 75-80 ns, while 
the nucleation of a new branch takes less than 5 ns (Fig. 5). The 
physical origin of these time scales can be understood by considering 
a combination of several physical processes, discussed one by one 
below. 

Firstly, the adiabatic heating time x a describes the shortest possible 
time to increase the temperature from T 0 to T c . This characteristic 
time is estimated by considering Ohmic Joule heating, and solving 
cT = p 0 j 2 c (T 0 ), where c = c(T c )(T/T c ) 3 is the specific heat. Integrating 
this equation gives 

r a = c(T c )T c /4 Po j 2 c (T 0 ), (1) 

where a small term ~ (T 0 /T c ) 4 is ignored. Using the material para- 
meters given in the Methods section, the numerical value becomes x a 
= 0.5 ns, which indeed is very close to the rise time of the simulated 
events. 

Secondly, the electromagnetic time scale T em describes the life- 
time of normal currents. For a thin disk, Brandt has found that the 
longest surviving mode has a decay time given by 39 

Tem = 0.18^ 0 Rd/p Q . (2) 

With the present parameters, this gives T em = 1.9 ns. Note that in the 
bulk case such a time constant cannot be defined since the flux 
motion is then described by a diffusion equation. In films, on the 
other hand, the flux penetration is accelerated by the presence of a 
free surface. The decay time is related to the propagation velocity of 
the peak in the current density, which is v em = 0.77p o /in o d = 0.14R/ 
i em = 140 km/s 39 . This value provides the upper bound for the 
propagation velocity of the dendrite. Indeed, the initial dendrite tip 
velocity ~ 90 km/s of Fig. 6 is not far from v em . 

Thirdly, heat removal to the substrate leads to an exponential 
decay of the temperature with a time constant 

Th = cd/h = 52 ns, (3) 

where h is the coefficient of heat transfer to the substrate. This decay 
time is much longer than both i a and T em , fully consistent with the 
fact that the events actually do take place, rather than being prevented 
by an efficient heat sink provided by the substrate. The value of x h is 
comparable to the total duration of the event, suggesting that the heat 
removal to the substrate largely determines the avalanche life-time, 
and thereby also decides the size of the full-grown flux dendrite. 

Finally, the lateral heat transport is an ordinary diffusion process 
with diffusion time 

t k = 1 2 c/4k, (4) 

where / is the diffusion length and k is the thermal conductivity. The 
diffusion length characteristic for the dendrite tips can be obtained 



from the T-maps of Fig. 5, where the very sharp tips of the growing 
branches have a typical width / = 10 /mi, which gives t k = 5.1 ns. 
This is close to the characteristic time for the first branching to occur, 
indicating that the heat diffusion plays an important role for the way 
the branch width evolves. 

Considering the other extreme, and letting t k be the total duration 
of an event, 75 ns, we obtain the largest relevant diffusion length, / = 
38 fim. This is much smaller than the length of the long branches in 
the dendritic structure, but interestingly it is approximately half the 
width of the dendrite trunk at the final stage. This indicates that the 
trunk is gradually widened by heat diffusion during the event. 

All these characteristic times are much shorter than the time scale 
of the background flux penetration, which is on the order of milli- 
seconds. Therefore, our results on the evolution of the instability are 
essentially independent of the ramp rate of the applied magnetic 
field. This robustness is consistent with numerous MOI experiments 
performed by some of the present authors. 

Comparison of the characteristic times shows that the propagation 
of flux avalanches in superconducting films placed in a transverse 
applied magnetic field is at the initial stage driven by nonlocal elec- 
trodynamics and adiabatic heating rather than by thermal diffusion. 
The thermal aspects of the present phenomenon have similarities to a 
moving hot spot generated by a laser beam swept across a surface at 
supersonic speed. The temperature profile is then determined by 
adiabatic heating from the fast moving source. Very fast flux propaga- 
tion due to adiabatic heating by the electromagnetic field outside a 
superconducting films was first pointed out in Ref. 28. Since x a <C t k , 
in the initial (supersonic) stage of flux propagation the thermal dif- 
fusion plays only a minor role. During later (subsonic) stages, thermal 
diffusion in the avalanche propagation direction will only slightly 
smear the temperature distribution formed by the local adiabatic 
heating. In the direction perpendicular to the avalanche, however, 
the thermal diffusion is indeed important as it leads to branching 
and determines the width of the branches. We have confirmed this 
by performing additional simulations for different values of k. 

In conclusion, by using a numerical method we have simulated the 
coupled dynamics of magnetic flux, electrical field and local temper- 
ature in superconducting films placed in a time- dependent trans- 
verse applied magnetic field. The calculations were performed for a 
disk-shaped thin sample with a set of parameters representative for 
MgB 2 , a material where breakdown events are frequently observed, 
but so far no experiments could resolve the ultra-fast dynamics. 

The temporal characteristics of the events are extreme, as the flux 
front propagates with the speed of 100 km/s and the whole event is 
over in less than 100 ns. We have found that the dendritic features 
first appear approximately 5 ns after the instability onset, when the 
root splits into branches that propagate into the flux-free interior of 
the sample. From that moment on, the morphology of the flux pat- 
tern is essentially self- similar, and maintained by side branches 
growing from the main branches as the structure expands. Occa- 
sionally, the branch tips will bifurcate symmetrically. The maximum 
temperature is found at the main trunk of the dendritic structure, 
where it reaches (1.2 - 1.4)T C , obviously fatal for the local super- 
conducting properties, and extremely harmful for a device based 
on such a superconducting film. In the same region, also the electric 
field is maximum, with values close to the theoretical maximum, 
Pojc(To) = 6.7 kV/m, reached almost immediately, f<l ns, after 
nucleation of the event. The maximum dissipation rate is approxi- 
mately 1 W. 

Prior to an event we find large fluctuations in both temperature 
and electric field. They are interpreted as traces of micro- avalanches, 
a phenomenon that so far has attracted little attention. We argue that 
these fluctuations are of prime importance, since they may well be the 
stochastic triggering mechanism causing the large-scale breakdown. 

Our results are in good agreement with the very few experimental 
works reporting on the time-resolved intermittent behavior in 
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superconducting films. Since adequate experimental tools are almost 
non-existent, our simulations method can be essential for the further 
exploration of the nucleation and evolution of these dramatic irre- 
producible and unpredictable lightning-like events. 

Methods 

Basic equations and numerical approach. The sample is a disk with radius R and 
thickness d <C placed in an applied magnetic field H a directed transverse to the disk 
plane, see Fig. 1 . We treat the magnetic flux density B z as a continuum variable 
satisfying the Maxwell equations, and assume a vanishing first critical field, H cX = 0, 
implying that B = ^ 0 H also inside the superconductor 35 . 

The material characteristics is modelled by the commonly used phenomenological 
non-retarded E - J relation 25 , 

£= PJ p fPo J>L(T), 
d ' XpoU/JcT' 1 otherwise. 

Here E is the electric field, J is the sheet current with |J| = /, p 0 is the normal state 
resistivity, J c is the critical sheet current, and n is the creep exponent. For T < T c the 
temperature dependencies of J c and n are specified as 

h=Jco(l-T/T c ), n = n 1 (T c /T)-n 0 . (6) 

This model describes thermally activated hopping of vortices, the effective velocity 
being oc e~ u(r)lkT with U = U 0 ln(/ c //). The model captures main features of nonlinear 
vortex dynamics and turns out to be very successful in understanding results of 
numerous experiments, see, e.g., Ref. 41 for a review. 

The sheet current density is expressed through the local magnetization, g(x, y), as 25 

dg/dy=J x , dg/dx=-J y , (7) 

which ensures conservation of current, V • J = 0. As described in Ref. 32, the time 
derivative of g can be obtained by inversion of the Biot-Savart law, which is expressed 
by 

g = T- l [(2/k)T{B z /n 0 -H a )]. (8) 

Here T and T~ 1 denote, respectively, the direct and inverse Fourier transforms with 
respect to spatial coordinates, and k = |k| where k is the in-plane wave vector. In Eq. 

(8) the three-dimensional Maxwell equations have been transformed to a non-local 
two-dimensional equation. The non-locality is a consequence of the in-plane mag- 
netic fields present at both surfaces of the film. The non-local electrodynamics of a 
thin film in the transverse geometry is the reason for fast propagation of flux dendrites 
with characteristic speeds exceeding microscopic electron and phonon scales such as 
the sound velocity. 

The time evolution of the temperature is modelled by the heat propagation 
equation, 

cf = V-(icVT) - (h/d)(T- T 0 ) + VE/d. (9) 

Here T 0 is the substrate temperature assumed to be constant, c is the specific heat. In 
the simulations we assume that c, h and k are all proportional to T 3 . These depen- 
dences are compatible with available experimental results 30 ' 40 . Since the characteristic 
phonon times are much shorter than the typical time scales for the development of the 
avalanche we neglect the delay in the thermal response, i. e., assume that both k and h 
do not explicitly depend on time. 

Equation (9) assumes that both electrons and phonons are characterized by quasi- 
equilibrium distributions with the same temperature. This is expected when the 
microscopic relaxation times corresponding to electron -electron, electron-phonon 
and phonon-phonon scattering are much shorter than the typical evolution times of 
the present thermo- magnetic propagation process. In addition, the mean free paths of 
the electrons and phonons should be much less than typical dimensions of the 
dendritic patterns. Estimates show that these requirements are met in the present 
situation. 

The evolution of g and T is computed by discrete time integration of Eqs. (8) and 

(9) . The plane z = 0 is arranged as a square superlattice where each cell has size 2a X 
2a, with a> R. The superconducting disk is placed centrally inside the cell and 
surrounded by vacuum. Then, within the superconductor, B z is found from Faraday's 
law, B z =-(Vx E) z , combined with the material law, Eq. (5). In the surrounding 
vacuum, B z is obtained by a self- consistent iterative scheme which ensures 

g = 0 there. Eq. (8) was integrated forward in time with the fourth order Runge-Kutta 
method, while the Eq. (9) was integrated by a forward/backward time average to 
ensure stable results. A detailed description of the computational scheme is found in 
Ref. 32. 

The simulations were run on a square of half- width a = 1.2R, and discretized on a 
512 X 512 equidistant grid. A weak quenched disorder was introduced by randomly 
increasing or decreasing / c0 of each grid point by 5%. In addition, fluctuations were 
introduced by randomly alternating between the left and the right derivatives when 
computing V X E in the Faraday law. 

Material parameters. The following parameters, typical for films of MgB 2 , were used 
in the calculations 30 ' 32 ' 40 : T c = 39 K, p 0 = 7 fiQ'cm, and a critical current density j c0 = 



] c Jd = 1.2 • 10 11 Aim 2 . The thermal parameters were chosen as k = 0.17 (T/T c ) 3 kW/ 
Km, c = 35 {TIT c f kJ/Km 3 , and h = 200 (T/T c ) 3 kW/Km 2 . We used n x = 20, n Q = 10, 
and restricted n(T) to n < 59 at low temperatures. The disk dimensions were R = 
2 mm and d = 0.3 fim. Substrate temperature was T 0 = 0.2T C and the applied 
magnetic field was ramped at the constant rate, H a = 2A-lO~ 6 j c op Q /Rfi 0 . 
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